install.packages("WRS", repos="http://R-Forge.R-project.org")
library(WRS)
source("http://www.statpower.net/Content/311/R Stuff/WRS.addon.txt")
js.t.test <- function(data1,data2,alpha){
  n1 <- length(data1)
  n2 <- length(data2)
  df <- n1 + n2 - 2
  crit.t <- qt(1-alpha/2,df)
  sigma.hat.squared <- ((n1-1)*var(data1) + (n2-1)*var(data2))/df
  t.obtained <- (mean(data1)-mean(data2))/sqrt(((1/n1)+(1/n2))*sigma.hat.squared)
  return(c(t.obtained,abs(t.obtained) > crit.t))
}
set.seed(12345)
rowMeans(replicate(10000,js.t.test(rnorm(10,0,sqrt(40)),rnorm(20,0,sqrt(10)),.05)))

make.data <-function(){
  y <- c(rnorm(100,50,10),rnorm(20,50,30))
  group <- c(rep(1,100),rep(2,20))
  data <- cbind(y,group)
  return(data)
}
tdata<-make.data()
t.test(y~group,var.equal=TRUE,data=tdata)

t.test(y~group,data=tdata)

output <- t.test(y~group,var.equal=TRUE,data=make.data())
str(output)

set.seed(12345)
replicate(5,t.test(y~group,var.equal=TRUE,data=make.data())$p.value)
set.seed(12345)
replicate(5,t.test(y~group,var.equal=TRUE,data=make.data())$p.value<.05)

set.seed(12345)
mean(replicate(2000,t.test(y~group,
                           var.equal=TRUE,data=make.data())$p.value<.05))
set.seed(12345)
mean(replicate(2000,t.test(y~group,
                           var.equal=FALSE,data=make.data())$p.value<.05))

make.data <-function(n1,n2,mean1,mean2,sd1,sd2){
  y <- c(rnorm(n1,mean1,sd1),rnorm(n2,mean2,sd2))
  group <- c(rep(1,n1),rep(2,n2))
  data <- cbind(y,group)
  return(data)
}

rejection.rate <- function(n1,n2,mean1,mean2,sd1,sd2,
             reps=2000,nominal.alpha=0.05,var.equal=FALSE) {
  mean(replicate(reps,
       t.test(y~group,
       var.equal=var.equal,
       data=make.data(n1,n2,mean1,mean2,sd1,sd2)
       )$p.value<nominal.alpha))
}

set.seed(12345)
rejection.rate(100,20,50,50,10,30)
rejection.rate(45,5,0,0,5,50)